What the code does:

  • Calculates all stats for song sound pressure level (i.e. amplitude) covariates

Next steps

  • check if magnitude of response to treatments in significant models is explained by body size of lifting power

Load packages

## vector with package names
x <- c( "pbapply", "parallel", "ggplot2", "warbleR", "Rraven", "viridis", "readxl", "rptR", "MCMCglmm", "MuMIn", "corrplot", "lme4", "grid", "gridExtra")

aa <- lapply(x, function(y) {
  
  # check if installed, if not then install 
  if (!y %in% installed.packages()[,"Package"]) 
    install.packages(y) 

  # load package
  try(require(y, character.only = T), silent = T)
})

Functions and parameters

#functions and parameters
knitr::opts_knit$set(root.dir = normalizePath(".."))

knitr::opts_chunk$set(dpi = 50, fig.width = 12, warning = FALSE, message = FALSE) 

# ggplot2 theme
# theme_set(theme_classic(base_size = 20))

cut_path <- "./data/raw/cuts"

treatments <- c("Calibration", "Regular_singing", "Coordination", "After_chase",  
 "Before_playback", "After_playback", "Before_interaction", "After_interaction", "Before_noise", "After_noise")

# iterations for MCMCglmm models
itrns <- 100000

# functions from https://rdrr.io/rforge/rptR/src/R/rpt.mcmcLMM.R
rpt.mcmcLMM <- function(y, groups, CI=0.95, prior=NULL, verbose=FALSE, ...){
    # initial checks
    if(length(y)!= length(groups)) stop("y and group are of unequal length")
    # preparation
    groups <- factor(groups)
    if(is.null(prior)) prior <- list(R=list(V=1,n=10e-2), G=list(G1=list(V=1,n=10e-2)) )
    # point estimation according to model 8 and equation 9
    mod   <- MCMCglmm(y ~ 1, random=~groups, family="gaussian", data=data.frame(y=y,groups=groups), prior=prior, verbose=verbose, ...)
    var.a <- mod$VCV[,"groups"]
    var.e <- mod$VCV[,"units"]
    postR <- var.a / (var.a + var.e)
    # point estimate
    R     <- posterior.mode( postR )
    # credibility interval estimation from paterior distribution
    CI.R    <- coda::HPDinterval(postR,CI)[1,]
    se      <- sd(postR)
    # 'significance test'
    P     <- NA
    res = list(call=match.call(), datatype="Gaussian", method="LMM.MCMC", CI=CI, 
                R=R, CI.R=CI.R, se=se, P=P) 
    # class(res) <- "rpt"
    return(res) 
}


## print Gelman-Rubin convergence statistics, plots traces and autocorrelations
mcmc_diagnostics <- function(rep_mods_list){
  
  for(w in 1:length(rep_mods_list))
{
  mod_name <- names(rep_mods_list)[w]
  
   if(mod_name == "1") mod_name <- "Null" 
  print(paste("model:", mod_name))
  
  Y <- lapply(rep_mods_list[[w]], "[[", "Sol")
  

  ## add global plots and gelman test
  # gelman_diagnostic
  gel_diag <- as.data.frame(gelman.diag(mcmc.list(Y))$psrf)
  
  # add estimate as column
  gel_diag$estimate <- rownames(gel_diag)
  
  # reorder columns
  gel_diag <- gel_diag[, c(3, 1, 2)]

  par(mfrow = c(1, 1))

  # plot table
  grid.newpage()
  grid.draw(tableGrob(gel_diag, rows = NULL, theme=ttheme_default(base_size = 25)))  

  par(mfrow = c(1, 4))

  traceplot(Y, col = adjustcolor(c("yellow","blue", "red"), alpha.f = 0.6))
  
  autocorr.plot(x = Y[[1]], auto.layout = FALSE, lwd =4, col = "red")
} 
  
}
amp <- read.csv("./output/calibrated_amplitude_all_songs.csv")
amp$Treatment[amp$Treatment == "Regular_sining"] <- "Regular_singing"

Repeatability

We estimated repeatability using linear mixed models with MCMC sampling (MCMCglmm() function)

# keep only morning recs, only regular singing and only from 2021
amp_morn <- amp[amp$period == "morning" & amp$year == 2021 & amp$Treatment == "Regular_singing", ]

# keep only males recorded at least twice 
count_sf <- table(amp_morn$ID[!duplicated(amp_morn$org.sound.file)])

amp_morn_rep <- amp_morn[amp_morn$ID %in% names(count_sf)[count_sf > 1], ]

max_quantile <- c(0, 1, seq(10, 90, by = 10)) / 100
outliers <- c(0, 1, 2, 5, 10, 20) / 100

rep_grid <- expand.grid(max_quantile = max_quantile, outliers = outliers, only.low.outliers = c(TRUE, FALSE))

repts <- pblapply(1:nrow(rep_grid), cl = 6, function(x){  
  
  quant_subet <- lapply(unique(amp_morn_rep$org.sound.file), function(y){
    X <- amp_morn_rep[amp_morn_rep$org.sound.file == y, ]
    
    # remove outliers
    outlier_quant <- quantile(X$cal.spl, c(rep_grid$outliers[x], 1- rep_grid$outliers[x]))
  
        if (!rep_grid$only.low.outliers[x])
    X <- X[X$cal.spl >= outlier_quant[1] & X$cal.spl <= outlier_quant[2],] else  X <- X[X$cal.spl >= outlier_quant[1],]
    
    # quantlie for each max quantile
    quant <- quantile(X$cal.spl, probs = 1 - rep_grid$max_quantile[x])
    
    # subset
   X <- X[X$cal.spl >= quant, ]
  })
  
  quant_subet <- do.call(rbind, quant_subet)
  
  # frequentist
    # out <- rpt(cal.spl ~ (1 | ID), grname = "ID", data = quant_subet, datatype = "Gaussian", npermut = 0, nboot = 100)

  # bayesian
    out <- rpt.mcmcLMM(y = quant_subet$cal.spl, groups = quant_subet$ID, nitt = itrns)
    
  out_df <- data.frame(rep_grid$max_quantile[x], rep_grid$outliers[x], rep_grid$only.low.outliers[x], out$R, out$CI.R[1], out$CI.R[2])
  
  return(out_df)
  })

repts_df <- do.call(rbind, repts)

names(repts_df) <- c("max_quantile", "outliers", "only.low.outliers", "repeatability", "lowCI", "hiCI")

write.csv(repts_df, "./output/repeatability_optimization.csv", row.names = FALSE)
  • Different subsets of data for upper quantiles (x axis)
  • Pannels show the proportion of outliers removed
repts_df <- read.csv("./output/repeatability_optimization.csv")

pd <- position_dodge(width = 0.1)
ggplot(data = repts_df, aes(x = 1 - max_quantile, y = repeatability, color = only.low.outliers, group = only.low.outliers)) + 
  geom_hline(yintercept = 0.5, col = adjustcolor("red", alpha.f = 0.5)) +
  geom_point(size = 2, position = pd) +
  geom_errorbar(width=.05, aes(ymin = lowCI, ymax = hiCI), position = pd) +
  scale_color_viridis(discrete = TRUE, begin = 0.2, end = 0.8, alpha = 0.7) +
  geom_line(position = pd) + 
  labs(y = "Repeatability",  x = "Upper quantile used") +
  ylim(c(0, 1)) + xlim(c(1, 0)) +
  facet_wrap(~ outliers, scales = "fixed") + 
  theme_classic(base_size = 24)

repts_df$range <- repts_df$hiCI - repts_df$lowCI

# repts_df[order(repts_df$repeatability), c("max_quantile", "outliers", "only.low.outliers", "repeatability", "range")]

repts_df <- repts_df[order(repts_df$repeatability, decreasing = TRUE), ]

kable(head(repts_df))
max_quantile outliers only.low.outliers repeatability lowCI hiCI range
126 0.3 0.20 FALSE 0.6159727 0.2920327 0.8264643 0.5344317
94 0.4 0.02 FALSE 0.5707252 0.2818027 0.8171755 0.5353727
40 0.5 0.05 TRUE 0.5645533 0.2930169 0.8232646 0.5302477
98 0.8 0.02 FALSE 0.5600182 0.2848066 0.8065382 0.5217316
93 0.3 0.02 FALSE 0.5589800 0.2987828 0.8229877 0.5242049
127 0.4 0.20 FALSE 0.5588000 0.2928387 0.8197835 0.5269448
  • Excluding lowest values increases repeatability
  • Excluding outliers doesn’t have a strong effect but removing 0.1 outliers produced the highest repeatability
  • Removing lower tail outliers or both works similarly (we excluded only low tail outliers to keep more data)

Subsetting the 0.4 upper quartile after excluding 0.2 lower tail outliers:

# compose variable to remove low values and outliers based on repeatabiliy
amp$osf.treat <- paste(amp$org.sound.file, amp$Treatment, sep = "-")

rm_outlier_amp_l <- lapply(unique(amp$osf.treat), function(y){
    X <- amp[amp$osf.treat == y, ]
    
    # remove outliers
    outlier_quant <- quantile(X$cal.spl, c(repts_df$outliers[1], 1))
    X <- X[X$cal.spl >= outlier_quant[1] & X$cal.spl <= outlier_quant[2],]
    
    # quantlie for each max quantile (0.6 was selected due to high repeatability)
    quant <- quantile(X$cal.spl, probs = repts_df$max_quantile[1])
    
    # subset
   X <- X[X$cal.spl >= quant, ]
  })

rm_outlier_amp <- do.call(rbind, rm_outlier_amp_l)

Hypothesis testing analysis

We use bayesian MCMC generalized linear models with individual ID and song type as random factors (random intercepts) except indicated. An intercept-only (null) model is also included in the analyses. A loop is used to run these models. Each model is replicated three times with starting values sampled from a Z-distribution (“start” argument in MCMCglmm()). Diagnostic plots for MCMC model performance are shown at the end of this report.

SPL vs day period (morning / afternoon)

  • Only for individuals recorded in both periods
  1. Period of the day as predictor: \[SPL \sim + period + (1 | individual) + (1 | song\ type)\]
  2. Null model with no predictor (intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]
# subset data
period_data <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" & rm_outlier_amp$year == 2021 & rm_outlier_amp$ID %in% unique(rm_outlier_amp$ID[rm_outlier_amp$period == "afternoon"]), ]

# model SPL by body size
period_formulas <- c("cal.spl ~ 1", "cal.spl ~ period")

period_mods <- pblapply(period_formulas, function(x){
  
  replicate(n = 3, MCMCglmm(as.formula(x), random = ~ ID + songtype, data = period_data, verbose = FALSE, nitt = itrns,  start = list(QUASI = FALSE)), simplify = FALSE)
  
})

names(period_mods) <- gsub("cal.spl ~ ", "", period_formulas)
saveRDS(list(period_formulas = period_formulas, period_mods = period_mods), "./output/mcmcglmm_period_models.RDS")

Model selection

attach(readRDS("./output/mcmcglmm_period_models.RDS"))

mod_list <- lapply(period_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl ~ ", "", period_formulas)

model_selection <- model.sel(mod_list, rank="DIC")

model_selection
(Intercept) period family df logLik DIC delta weight
period 98.38209 + (NA) 5 -6192.009 12389.99 0.00000 1.000000e+00
1 99.28286 NA (NA) 4 -6223.692 12452.40 62.40826 2.806839e-14

Best model summary

# fixed effects with HPD intervals
best_mod_period <- period_mods[[which(names(period_mods) == row.names(model_selection)[1])]][[1]]

sm <- as.data.frame(summary(best_mod_period)$solutions[, -5])

sm
post.mean l-95% CI u-95% CI eff.samp
(Intercept) 98.382089 90.2424713 106.17904 4700
periodmorning 1.232936 0.9223919 1.52236 4700
dat_period <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" & rm_outlier_amp$ID %in% unique(rm_outlier_amp$ID[rm_outlier_amp$period == "afternoon"]), ]

dat_period$period <- factor(dat_period$period, levels = c("morning", "afternoon"))

ggplot(dat_period, aes(x = as.factor(ID), y = cal.spl, color = period, fill = period)) +
  geom_violin(position =  pd) +
  labs(x = "Individual", y = "Sound pressure level (dB)", fill = "Period") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.6) +
  scale_fill_viridis_d(begin = 0.1, end = 0.8, alpha = 0.6) +
  guides(color = FALSE) +
  theme_classic(base_size = 24) 

Effect size posterior distribution

# simplify names
colnames(best_mod_period$Sol) <- c("Intercept", "Morning vs afternoon")

# stack posteriors
Y <- stack(as.data.frame(best_mod_period$Sol[, -1]))
Y$ind <- "Morning vs afternoon"

# plot posteriors
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution",x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 0.9, hjust=1)) 

  • Morning songs have significantly higher SPL than afternoon songs within individuals

SPL vs condition

  • Only for individuals recorded in the morning (afternoon recordings were excluded)
  • Condition parameters: body size (PC1 from a PCA on total culmen, flattened wing length, body mass and central rectrix) and lifting power
  1. Body size as predictor: \[SPL \sim + body\ size + (1 | individual) + (1 | song\ type)\]

  2. Lifting power as predictor: \[SPL \sim + lifting\ power + (1 | individual) + (1 | song\ type)\]

  3. Lifting power and body size as predictors: \[SPL \sim + lifting\ power + body\ size + (1 | individual) + (1 | song\ type)\]

  4. Interaction of lifting power and body size as predictor: \[SPL \sim + lifting\ power * body\ size + (1 | individual) + (1 | song\ type)\]

  5. Null model with no predictor (intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

caps <- read_excel("/home/m/Dropbox/LBH data/Additional data files/LBH captures data.xlsx")

# same variables as in paper on spatial memory and body size  
size_vars <- caps[caps$`Bird ID` %in% amp$ID, c("Bird ID", "Total culmen", "Flattened wing length", "Weight", "Central rectriz")]

# replace 419 central rectriz with mean
size_vars$`Central rectriz`[size_vars$`Bird ID` == 419] <- mean(size_vars$`Central rectriz`, na.rm = TRUE)

cor(size_vars[, c("Total culmen", "Flattened wing length", "Weight", "Central rectriz")], use = "pairwise.complete.obs")
##                       Total culmen Flattened wing length     Weight
## Total culmen           1.000000000           -0.13969106 0.14426785
## Flattened wing length -0.139691059            1.00000000 0.08646515
## Weight                 0.144267849            0.08646515 1.00000000
## Central rectriz        0.007536892            0.34107800 0.10235729
##                       Central rectriz
## Total culmen              0.007536892
## Flattened wing length     0.341077996
## Weight                    0.102357292
## Central rectriz           1.000000000
mean_size <- aggregate(. ~ `Bird ID`, data = size_vars, FUN = mean, na.action = "na.omit")

pca <- prcomp(mean_size[, c("Total culmen", "Flattened wing length", "Weight", "Central rectriz")], scale. = TRUE)

summary(pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4
## Standard deviation     1.2488 1.0841 0.8726 0.7098
## Proportion of Variance 0.3899 0.2938 0.1903 0.1259
## Cumulative Proportion  0.3899 0.6837 0.8741 1.0000
mean_size$PC1 <- pca$x[, 1]


rm_outlier_amp$PC1 <- sapply(1:nrow(rm_outlier_amp), function(x) mean_size$PC1[mean_size$`Bird ID` == rm_outlier_amp$ID[x]][1])

rm_outlier_amp$weight <- sapply(1:nrow(rm_outlier_amp), function(x) mean_size$Weight[mean_size$`Bird ID` == rm_outlier_amp$ID[x]][1])

### add weight lifted
lift <- read_excel("/home/m/Dropbox/LBH data/Morphology-condition/other files/Uplift power experiment.xlsx", sheet = "Results")

rm_outlier_amp$lift.weight <- sapply(1:nrow(rm_outlier_amp), function(x){
  
  if(rm_outlier_amp$ID[x] == 0) weight <- NA else {
    
    if (rm_outlier_amp$year[x] == 2019)
      sub_lift <- lift[grep("\\.2019\\.", lift$Video), ] else
              sub_lift <- lift[grep("\\.2021\\.", lift$Video), ]

    
    weights <- sub_lift$Weight[sub_lift$ID == rm_outlier_amp$ID[x]]
   
    if (length(weights) > 0)
    # get the mean of the 2 highest flights
     weight <- mean(sort(weights, decreasing = TRUE)[1:2]) else
    weight <- NA
  }  
  
  return(weight)
})
# subset data
body_size_data <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" & rm_outlier_amp$year == 2021 & !is.na(rm_outlier_amp$PC1) & !is.na(rm_outlier_amp$lift.weight) & rm_outlier_amp$period == "morning", ]

# model SPL by body condition
size_formulas <- c("cal.spl ~ 1", "cal.spl ~ PC1", "cal.spl ~ lift.weight", "cal.spl ~ lift.weight + PC1", "cal.spl ~ lift.weight *  PC1")


size_mods <- pblapply(size_formulas, function(x){
  
  replicate(n = 3, MCMCglmm(as.formula(x), random = ~ ID + songtype, data = body_size_data, verbose = FALSE, nitt = itrns,  start = list(QUASI = FALSE)), simplify = FALSE)
  
})

names(size_mods) <- gsub("cal.spl ~ ", "", size_formulas)
saveRDS(list(size_formulas = size_formulas, size_mods = size_mods), "./output/mcmcglmm_size_models.RDS")

Model selection

attach(readRDS("./output/mcmcglmm_size_models.RDS"))

mod_list <- lapply(size_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl ~ ", "", size_formulas)

model_selection <- model.sel(mod_list, rank="DIC")

model_selection
(Intercept) PC1 lift.weight lift.weight:PC1 family df logLik DIC delta weight
PC1 101.06722 1.0132676 NA NA (NA) 5 -19465.51 38951.92 0.0000000 0.2082655
lift.weight + PC1 85.44995 0.8044757 0.9775365 NA (NA) 6 -19465.54 38951.92 0.0073858 0.2074978
lift.weight * PC1 77.28863 10.1490394 1.5058787 -0.6219876 (NA) 7 -19465.49 38952.02 0.1043244 0.1976803
lift.weight 84.06902 NA 1.0776599 NA (NA) 5 -19465.55 38952.06 0.1464912 0.1935562
1 101.35748 NA NA NA (NA) 4 -19465.57 38952.07 0.1522446 0.1930002

Best model summary

None of the models provided a better fit compared to the null model

# fixed effects with HPD intervals
best_mod_size <- size_mods[[which(names(size_mods) ==  "lift.weight + PC1")]][[1]]

sm_size <- as.data.frame(summary(best_mod_size)$solutions[, -5])

sm_size
post.mean l-95% CI u-95% CI eff.samp
(Intercept) 85.4499479 49.409360 122.753355 4700.00
lift.weight 0.9775365 -1.421323 3.169801 4700.00
PC1 0.8044757 -2.093328 3.909391 4354.76
body_size_data <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" & rm_outlier_amp$year == 2021 & !is.na(rm_outlier_amp$PC1) & !is.na(rm_outlier_amp$lift.weight) & rm_outlier_amp$period == "morning", ]

agg_size <- aggregate(cbind(cal.spl, lift.weight, PC1) ~ ID, data = body_size_data, mean)
agg_size$count <- aggregate(PC1 ~ ID, data = body_size_data, length)$PC1

ggplot(agg_size, aes(x = PC1, y = cal.spl)) +
    geom_jitter(position = position_jitter(width = 0.2, height = 0.2), aes(size = count)) +
  labs(x = "Body size (PC1)", y = "Sound pressure level (dB)", color = "Period", size = "Sample size") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24) + 
  # geom_abline(slope = c(ES_continious = sm_size[2, 1], ES_continious_plus_ES_interaction=  sm_size[2, 1] + sm_size[4, 1]), intercept = c(intercept = sm_size[1, 1], intercept_plus_ES_categorical = sm_size[1, 1] + sm_size[3, 1]), color = viridis(2, begin = 0.2, end = 0.8, alpha = 0.6)[2:1], size = 1) +
guides(colour = guide_legend(override.aes = list(size=8)))

ggplot(agg_size, aes(x = PC1, y = cal.spl)) +
    geom_jitter(position = position_jitter(width = 0.2, height = 0.2), aes(size = count)) +
  labs(x = "Body size (PC1)", y = "Sound pressure level (dB)", color = "Period", size = "Sample size") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24) + 
  # geom_abline(slope = c(ES_continious = sm_size[2, 1], ES_continious_plus_ES_interaction=  sm_size[2, 1] + sm_size[4, 1]), intercept = c(intercept = sm_size[1, 1], intercept_plus_ES_categorical = sm_size[1, 1] + sm_size[3, 1]), color = viridis(2, begin = 0.2, end = 0.8, alpha = 0.6)[2:1], size = 1) +
guides(colour = guide_legend(override.aes = list(size=8)))

ggplot(agg_size, aes(x = lift.weight, y = cal.spl)) +
    geom_jitter(position = position_jitter(width = 0.2, height = 0.2), aes(size = count)) +
  labs(x = "Lifiting power (g)", y = "Sound pressure level (dB)", color = "Period", size = "Sample size") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24) + 
  # geom_abline(slope = c(ES_continious = sm_size[2, 1], ES_continious_plus_ES_interaction=  sm_size[2, 1] + sm_size[4, 1]), intercept = c(intercept = sm_size[1, 1], intercept_plus_ES_categorical = sm_size[1, 1] + sm_size[3, 1]), color = viridis(2, begin = 0.2, end = 0.8, alpha = 0.6)[2:1], size = 1) +
guides(colour = guide_legend(override.aes = list(size=8)))

  • Each dot represents the average for a single individual

Effect size posterior distribution

# simplify names
colnames(best_mod_size$Sol) <- c("Intercept", "lifted weight","Body size (PC1)")

# stack posteriors
Y <- stack(as.data.frame(best_mod_size$Sol[, -1]))

# plot posteriors
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution", x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 1, hjust=1)) 

  • Song PSL is not affected by condition (either body size or lifting power)

SPL vs playback

  • Using data from 2019 and 2021
  • Recordings from 2019 were not calibrated and were center around SPL mean from 2021
  • Comparing before playback vs after playback (during regular singing)
  1. Playback as predictor and random slope: \[SPL \sim + body\ size + (playback | individual) + (1 | song\ type)\]

  2. Playback as predictor and random slope: \[SPL \sim + body\ size + (playback | individual) + (1 | song\ type)\]

  3. Playback as predictor and only random intercept: \[SPL \sim + playback + (1 | individual) + (1 | song\ type)\]

  4. Null model with no predictor (random intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

playback_dat <- rm_outlier_amp[rm_outlier_amp$Treatment %in% c("Before_playback", "After_playback"), ]

# label 2019 recordings
playback_dat$ID <- ifelse(grepl("2019", playback_dat$org.sound.file), paste(playback_dat$ID, "(2019)"), playback_dat$ID)

# remove indiviudal recorded in 2019 and 2021
playback_dat <- playback_dat[playback_dat$ID != "397 (2019)", ]

random_slope_mod <- lmer(cal.spl ~ Treatment  + (Treatment|ID), data = playback_dat, REML = FALSE)

interaction_size_mod <- lmer(cal.spl ~ PC1 + Treatment + PC1:Treatment + (1|ID), data = playback_dat, REML = FALSE)

interaction_lift_mod <- lmer(cal.spl ~ PC1 + lift.weight + PC1:lift.weight + (1|ID), data = playback_dat, REML = FALSE)

random_intercept_only_mod <- lmer(cal.spl ~ Treatment  + (1|ID) + (1|songtype), data = playback_dat, REML = FALSE)

null_mod <- lmer(cal.spl ~ 1  + (1|ID) + (1|songtype), data = playback_dat, REML = FALSE)
d1 <- (-2*logLik(random_intercept_only_mod)) + (2*logLik(random_slope_mod))

pval_rand_slope_vs_rand_intercept <- 0.5*pchisq(d1[1], 0, lower.tail=FALSE) + 0.5*pchisq(d1[1], 1, lower.tail=FALSE)

d2 <- (-2*logLik(null_mod)) + (2*logLik(random_intercept_only_mod))

pval_rand_intercept_vs_null <- 0.5*pchisq(d2[1], 0, lower.tail=FALSE) + 0.5*pchisq(d2[1], 1, lower.tail=FALSE)

d3 <- (-2*logLik(null_mod)) + (2*logLik(random_slope_mod))

pval_rand_slope_vs_null <- 0.5*pchisq(d3[1], 0, lower.tail=FALSE) + 0.5*pchisq(d3[1], 1, lower.tail=FALSE)

d4 <- (-2*logLik(interaction_size_mod)) + (2*logLik(random_slope_mod))

pval_rand_slope_vs_interaction_size <- 0.5*pchisq(d4[1], 0, lower.tail=FALSE) + 0.5*pchisq(d4[1], 1, lower.tail=FALSE)

d5 <- (-2*logLik(interaction_lift_mod)) + (2*logLik(random_slope_mod))

pval_rand_slope_vs_interaction_lift <- 0.5*pchisq(d5[1], 0, lower.tail=FALSE) + 0.5*pchisq(d5[1], 1, lower.tail=FALSE)
agg_playback <- aggregate(cal.spl ~ Treatment + ID, data = playback_dat, mean)

agg_playback$sd <- aggregate(cal.spl ~ Treatment + ID, data = playback_dat, sd)$cal.spl

agg_playback$Treatment <- factor(gsub("_playback", "", agg_playback$Treatment), levels = c("Before", "After"))

ggplot(agg_playback, aes(x = Treatment, y = cal.spl, color = Treatment)) +
  geom_point(size = 2, show.legend = FALSE) +
   geom_errorbar(width=.05, aes(ymin = cal.spl - sd, ymax = cal.spl + sd), show.legend = FALSE) +
  labs(x = "Playback treatment", y = "Sound pressure level (dB)", color = "Period") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 1) +
  facet_wrap(~ ID, nrow = 3, scales = "free_y") +
  theme_classic(base_size = 24)

  • The random slope model is significantly better than the random intercept model (p = < 0.0001), the interaction with body size model (p = < 0.0001), the interaction with lifting power model (p = < 0.0001) and the null model (p = < 0.0001)
  • The playback stimuli elicits a significant change in SPL but the direction of the response varies between individuals
  • Direction of the response doesn’t seem to be affected by body size

SPL vs coordination

  • Random slope models, using data from 2019 and 2021
  • Only for individuals in which coordination was observed when recorded
  • SPL of songs sang during coordinated singing were compared to regular singing songs
  1. Coordination as predictor and random slope: \[SPL \sim + coordination + (coordination | individual) + (1 | song\ type)\] 1.Coordination as predictor and random intercept only: \[SPL \sim + coordination + (1 | individual) + (1 | song\ type)\]

  2. Null model with no predictor (random intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

# fixed effects with HPD intervals
coord_data <- rm_outlier_amp[rm_outlier_amp$org.sound.file %in% rm_outlier_amp$org.sound.file[rm_outlier_amp$Treatment == "Coordination"], ]

coord_data <- coord_data[coord_data$Treatment %in%  c("Regular_singing","Coordination"), ]

random_slope_mod <- lmer(cal.spl ~ Treatment  + (Treatment|ID) + (1|songtype), data = coord_data, REML = FALSE)

random_intercept_only_mod <- lmer(cal.spl ~ Treatment + (1|ID) + (1|songtype), data = coord_data, REML = FALSE)

null_mod <- lmer(cal.spl ~ 1  + (1|ID) + (1|songtype), data = coord_data, REML = FALSE)
d1 <- (-2*logLik(random_intercept_only_mod)) + (2*logLik(random_slope_mod))

# 0.5*pchisq(d1[1], 0, lower.tail=FALSE) + 0.5*pchisq(d1[1], 1, lower.tail=FALSE)


d2 <- (-2*logLik(null_mod)) + (2*logLik(random_intercept_only_mod))

# 0.5*pchisq(d2[1], 0, lower.tail=FALSE) + 0.5*pchisq(d2[1], 1, lower.tail=FALSE)

d3 <- (-2*logLik(null_mod)) + (2*logLik(random_slope_mod))

# 0.5*pchisq(d3[1], 0, lower.tail=FALSE) + 0.5*pchisq(d3[1], 1, lower.tail=FALSE)

# sm_playback <- as.data.frame(fixef(random_slope_mod), confint(random_slope_mod))
# 
# sm_playback

summary(random_slope_mod)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: cal.spl ~ Treatment + (Treatment | ID) + (1 | songtype)
##    Data: coord_data
## 
##      AIC      BIC   logLik deviance df.resid 
##  21588.1  21631.8 -10787.1  21574.1     3768 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9918 -0.1442  0.0093  0.1648  3.8007 
## 
## Random effects:
##  Groups   Name                     Variance  Std.Dev. Corr 
##  ID       (Intercept)              7.513e+01 8.667822      
##           TreatmentRegular_singing 3.656e+00 1.912156 -0.26
##  songtype (Intercept)              4.264e-06 0.002065      
##  Residual                          1.730e+01 4.158752      
## Number of obs: 3775, groups:  ID, 11; songtype, 9
## 
## Fixed effects:
##                          Estimate Std. Error t value
## (Intercept)              101.7875     2.6173  38.890
## TreatmentRegular_singing  -0.2087     0.6293  -0.332
## 
## Correlation of Fixed Effects:
##             (Intr)
## TrtmntRglr_ -0.251
agg_coord <- aggregate(cal.spl ~ Treatment + ID, data = coord_data, mean)

agg_coord$sd <- aggregate(cal.spl ~ Treatment + ID, data = coord_data, sd)$cal.spl

agg_coord$Treatment <- factor(ifelse(agg_coord$Treatment == "Coordination", "Coordinated", "Solo"), levels = c("Solo", "Coordinated"))

ggplot(agg_coord, aes(x = Treatment, y = cal.spl, color = Treatment)) +
  geom_point(size = 2, show.legend = FALSE) +
   geom_errorbar(width=.05, aes(ymin = cal.spl - sd, ymax = cal.spl + sd), show.legend = FALSE) +
 
  labs(x = "coord treatment", y = "Sound pressure level (dB)", color = "Period") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 1) +
  facet_wrap(~ ID, nrow = 2, scales = "free_y") +
  theme_classic(base_size = 24)

  • The random slope model was not significantly better than the random intercept model (p = < 0.0001) nor the null model (p = < 0.0001)
  • Vocal coordination is not associated with detectable changes in SPL

SPL vs aggresive interaction

  • 2 types of interactions: perch interaction and chases
  • SPL of songs before (i.e. regular singing) and after interactions were compared
  • Interactions encoded in 2 ways:
      1. factor clumping perch interactions and chases in a single category (2 levels: before and after interaction)
      1. Interaction type (chase of perch interaction)
  1. Interaction as predictor: \[SPL \sim + interaction + (1 | individual) + (1 | song\ type)\]
  2. Interaction and interaction type as predictors: \[SPL \sim + interaction + interaction\ type + (1 | individual) + (1 | song\ type)\]
  3. Statistical interaction between “Interaction” and “interaction type” as predictor: \[SPL \sim + interaction * interaction\ type + (1 | individual) + (1 | song\ type)\]
  4. Null model with no predictor (intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]
# subset data
dat_chase <- rm_outlier_amp[rm_outlier_amp$org.sound.file %in% rm_outlier_amp$org.sound.file[rm_outlier_amp$Treatment == "After_chase"], ]

dat_chase <- dat_chase[dat_chase$Treatment %in% c("Regular_singing", "After_chase"),]


dat_chase$inter_treatment <- ifelse(grepl("After", dat_chase$Treatment), "After", "Before")

dat_chase$inter_type <- "Chase"

dat_interaction <- rm_outlier_amp[rm_outlier_amp$Treatment %in% c("After_interaction", "Before_interaction"), ]


dat_interaction$inter_treatment <- ifelse(grepl("After", dat_interaction$Treatment), "After", "Before")
dat_interaction$inter_type <- "Perch_interaction"

# fix range of non-calibrated recs
dat_interaction$cal.spl[dat_interaction$year < 2021] <- dat_interaction$cal.spl[dat_interaction$year < 2021] + mean(dat_interaction$cal.spl[dat_interaction$year == 2021]) - mean(dat_interaction$cal.spl[dat_interaction$year < 2021])

dat_aggresive_inter <- rbind(dat_interaction, dat_chase)

# fix low value one
dat_aggresive_inter$cal.spl[dat_aggresive_inter$sound.files == dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl)]] <- dat_aggresive_inter$cal.spl[dat_aggresive_inter$sound.files == dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl)]]  +  mean(dat_aggresive_inter$cal.spl[dat_aggresive_inter$sound.files != dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl)]]) - mean(dat_aggresive_inter$cal.spl[dat_aggresive_inter$sound.files == dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl)]])

# table(dat_aggresive_inter$ID, dat_aggresive_inter$inter_treatment, dat_aggresive_inter$inter_type)
# model SPL by interaction
interaction_formulas <- c("cal.spl ~ 1", "cal.spl ~ inter_treatment", "cal.spl ~ inter_treatment + inter_type", "cal.spl ~ inter_treatment * inter_type")

inter_mods <- pblapply(interaction_formulas, function(x){
  
  replicate(n = 3, MCMCglmm(as.formula(x), random = ~ ID + songtype, data = dat_aggresive_inter, verbose = FALSE, nitt = itrns,  start = list(QUASI = FALSE)), simplify = FALSE)
  
})

names(inter_mods) <- gsub("cal.spl ~ ", "", interaction_formulas)

saveRDS(list(interaction_formulas = interaction_formulas, inter_mods = inter_mods), "./output/mcmcglmm_interaction_models.RDS")

Model selection

attach(readRDS("./output/mcmcglmm_interaction_models.RDS"))

mod_list <- lapply(inter_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl ~ ", "", interaction_formulas)

model_selection <- model.sel(mod_list, rank="DIC")

model_selection
(Intercept) inter_treatment inter_type inter_treatment:inter_type family df logLik DIC delta weight
inter_treatment * inter_type 105.8312 + + + (NA) 7 -2926.655 5874.963 0.00000 1.000000e+00
inter_treatment + inter_type 104.1285 + + NA (NA) 6 -2960.277 5941.127 66.16386 4.292404e-15
inter_treatment 102.7296 + NA NA (NA) 5 -2973.693 5966.872 91.90835 1.102439e-20
1 101.4128 NA NA NA (NA) 4 -3051.179 6120.967 246.00364 3.810561e-54

Best model summary

# fixed effects with HPD intervals
best_mod_interaction <- inter_mods[[which(names(inter_mods) == row.names(model_selection)[1])]][[1]]

sm_interaction <- as.data.frame(summary(best_mod_interaction)$solutions[, -5])

sm_interaction
post.mean l-95% CI u-95% CI eff.samp
(Intercept) 105.831196 103.633244 108.344906 1993.982
inter_treatmentBefore -4.252822 -4.765444 -3.725284 4700.000
inter_typePerch_interaction -3.970494 -4.775029 -3.131514 4700.000
inter_treatmentBefore:inter_typePerch_interaction 3.042986 2.351168 3.764623 4700.000
# subset data
dat_aggresive_inter$inter_treatment <- factor(dat_aggresive_inter$inter_treatment, levels = c("Before", "After"))

dat_aggresive_inter$year <- factor(dat_aggresive_inter$year, levels = c("2013", "2014", "2019", "2021"))

dat_aggresive_inter$inter_type <- gsub("_", " ", dat_aggresive_inter$inter_type)
  
dat_aggresive_inter$inter_type <- factor(dat_aggresive_inter$inter_type, levels = c("Perch interaction", "Chase"))

agg_aggresive_inter <- aggregate(cal.spl ~ ID + inter_type + inter_treatment + year, data = dat_aggresive_inter, mean)

agg_aggresive_inter$ID <- as.character(agg_aggresive_inter$ID)

ggplot(dat_aggresive_inter, aes(x = inter_treatment, y = cal.spl)) +
  geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  geom_point(size = 4, data = agg_aggresive_inter, aes(group = ID, color = year)) +
  geom_line(data = agg_aggresive_inter, aes(group = ID, color = year)) +
  labs(x = "Context", y = "Sound pressure level (dB)") +
  facet_wrap(~ inter_type) +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24) 

Effect size posterior distribution

# simplify names
colnames(best_mod_interaction$Sol) <- c("Intercept", "After vs before interaction", "Chase vs perch_interaction", "After vs before chase")

# stack posteriors
Y <- stack(as.data.frame(best_mod_interaction$Sol[, -1]))

# plot posteriors
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution",x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 0.9, hjust=1)) 

  • Male-male aggressive interactions result in a significant SPL increase
  • Chases have a stronger effect, increasing SPL in ~ 4 dB compare to regular singing before the chase

SPL and song structure

Song structure measured as spectrographic consistency and gap duration

Song spectrographic consistency

  • spectrographic consistency: mean spectrographic cross-correlation of all songs in a singing bout
  • Spectrographic consistency was measured for all songs within contexts (before & after interaction, clumping chase and perch interactions)
  • Spectrographic consistency was compared between the contexts
  • Recordings from interactions were used as they showed the highest increase in SPL
  1. interaction as predictor and spectrographic consistency as response: \[consistency \sim + interaction + (1 | individual) + (1 | song\ type)\]
  2. interaction and SPL as predictors and spectrographic consistency as response: \[consistency \sim + interaction + SPL + (1 | individual) + (1 | song\ type)\]
  3. interaction between “interaction” and SPL as predictor and spectrographic consistency as response: \[consistency \sim + interaction * SPL + (1 | individual) + (1 | song\ type)\]
  4. Null model with no predictor (intercept only): \[consistency \sim + 1 + (1 | individual) + (1 | song\ type)\]
agg_amp_inter <- aggregate(cal.spl ~ ID + inter_type, data = dat_aggresive_inter[dat_aggresive_inter$inter_treatment == "After",], mean)  

agg_amp_inter$amp.change <- agg_amp_inter$cal.spl -  
aggregate(cal.spl ~ ID + inter_type, data = dat_aggresive_inter[dat_aggresive_inter$inter_treatment == "Before",], mean)$cal.spl  

# sort by change
agg_amp_inter <- agg_amp_inter[order(agg_amp_inter$amp.change, decreasing = TRUE), ]

# subset data
dat_song_consistency <- dat_aggresive_inter[dat_aggresive_inter$ID %in% agg_amp_inter$ID[agg_amp_inter$amp.change >= 0.5], ]
dat_song_consistency$ID.treatment <- paste(dat_song_consistency$ID, dat_song_consistency$inter_treatment, sep ="-")

dat_song_consistency_l <- lapply(unique(dat_song_consistency$ID.treatment), function(x){
  
  X <- dat_song_consistency[dat_song_consistency$ID.treatment == x, ]
  
  print(x)
  if(nrow(X) > 1){
        # measure cross-correlation
    xcrr <- (cross_correlation(X, path = "./data/raw/cuts", pb = FALSE, wl = 150))  
    
    # measure gaps
    gps <- gaps(X, pb = FALSE)$gaps
    gps <- gps[gps < 1]
    mean.gaps <- mean(gps, na.rm = TRUE)
    
    # mean xcorr
    if (!is.na(xcrr[1])){
    xcrr <- mean(xcrr[lower.tri(xcrr)])
  }
  
  out_df <- data.frame(sound.files = x, ID = X$ID[1], Treatment = X$Treatment[1], xcrr = xcrr,  mean.gaps = mean.gaps, mean.spl = mean(X$cal.spl), songtype = X$songtype[1], treatment = X$inter_treatment[1])
  } else 
  out_df <- data.frame(sound.files = x, ID = X$ID[1], Treatment = X$Treatment[1], xcrr = NA,  mean.gaps = NA, mean.spl = mean(X$cal.spl), songtype = X$songtype[1], treatment = X$inter_treatment[1])
  
  return(out_df)
  
})

dat_song_consistency <- do.call(rbind, dat_song_consistency_l)

dat_song_consistency <- dat_song_consistency[!is.na(dat_song_consistency$xcrr), ]

saveRDS(dat_song_consistency, "./output/dat_song_spectrographic_consistency_and_gaps.RDS")
dat_song_consistency <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")


# model SPL by body size
song_consistency_formulas <- c("xcrr ~ 1", "xcrr ~ treatment", "xcrr ~ treatment + mean.spl", "xcrr ~ treatment * mean.spl" , "xcrr ~mean.spl")

song_consistency_mods <- pblapply(song_consistency_formulas, function(x){
  
  replicate(n = 3, MCMCglmm(as.formula(x), random = ~ ID + songtype, data = dat_song_consistency, verbose = FALSE, nitt = itrns,  start = list(QUASI = FALSE)), simplify = FALSE)
  
})

names(song_consistency_mods) <- gsub("cal.spl ~ ", "", song_consistency_formulas)

saveRDS(list(song_consistency_formulas = song_consistency_formulas, song_consistency_mods = song_consistency_mods), "./output/mcmcglmm_song_consistency_models.RDS")

Model selection

attach(readRDS("./output/mcmcglmm_song_consistency_models.RDS"))

mod_list <- lapply(song_consistency_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl ~ ", "", song_consistency_formulas)

model_selection <- model.sel(mod_list, rank="DIC")

model_selection
(Intercept) treatment mean.spl mean.spl:treatment family df logLik DIC delta weight
xcrr ~ treatment 0.7383678 + NA NA (NA) 5 34.73809 -62.44242 0.000000 0.62783916
xcrr ~ treatment + mean.spl 0.2875213 + 0.0043802 NA (NA) 6 32.91965 -59.39476 3.047666 0.13679057
xcrr ~ 1 0.7546914 NA NA NA (NA) 4 32.56895 -59.33922 3.103208 0.13304401
xcrr ~mean.spl 0.5292037 NA 0.0022100 NA (NA) 5 31.75161 -57.88083 4.561592 0.06416694
xcrr ~ treatment * mean.spl 0.2022626 + 0.0051972 + (NA) 7 31.84366 -56.84140 5.601029 0.03815931

Best model summary

# summary
best_mod_song_consistency <- song_consistency_mods[[which(names(song_consistency_mods) == row.names(model_selection)[1])]][[1]]

sm_song_consistency <- as.data.frame(summary(best_mod_song_consistency)$solutions[, -5])

sm_song_consistency
post.mean l-95% CI u-95% CI eff.samp
(Intercept) 0.7383678 0.6867783 0.7959804 4689.774
treatmentBefore 0.0323433 -0.0122364 0.0807822 9700.000
dat_song_consistency <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")

# subset data
dat_song_consistency$ID <- as.character(dat_song_consistency$ID)

ggplot(dat_song_consistency, aes(x = mean.spl, y = xcrr)) +
  geom_point(size = 4, color = viridis(5, alpha = 0.5)[3]) +
  labs(y = "Mean cross-correlation", x = "Sound pressure level (dB)") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24) 

Effect size posterior distribution

# stack posteriors
Y <- stack(as.data.frame(best_mod_song_consistency$Sol[, -1]))
Y$ind <- "mean.spl"

# plot posteriors
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution",x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 0.9, hjust=1)) 

  • Aggressive interaction but not SPL affects the spectrographic consistency of songs

Song gap duration

Gaps: silent time intervals between songs

  1. interaction as predictor and gap duration as response: \[gap\ duration \sim + interaction + (1 | individual) + (1 | song\ type)\]
  2. interaction and SPL as predictors and gap duration as response: \[gap\ duration \sim + interaction + SPL + (1 | individual) + (1 | song\ type)\]
  3. interaction between “interaction” and SPL as predictor and gap duration as response: \[gap\ duration \sim + interaction * SPL + (1 | individual) + (1 | song\ type)\]
  4. Null model with no predictor (intercept only): \[gap\ duration \sim + 1 + (1 | individual) + (1 | song\ type)\]
dat_gaps <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")

# model SPL by body size
gap_formulas <- c("mean.gaps ~ 1", "mean.gaps ~ inter_treatment", "mean.gaps ~ inter_treatment + mean.spl", "mean.gaps ~ inter_treatment * mean.spl" , "mean.gaps ~ mean.spl")

gap_mods <- pblapply(gap_formulas, function(x){
  
  replicate(n = 3, MCMCglmm(as.formula(x), random = ~ ID + songtype, data = dat_gaps, verbose = FALSE, nitt = itrns,  start = list(QUASI = FALSE)), simplify = FALSE)
  
})

names(gap_mods) <- gsub("cal.spl ~ ", "", gap_formulas)

saveRDS(list(gap_formulas = gap_formulas, gap_mods = gap_mods), "./output/mcmcglmm_gap_models.RDS")

Model selection

attach(readRDS("./output/mcmcglmm_gap_models.RDS"))

mod_list <- lapply(gap_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl ~ ", "", gap_formulas)

model_selection <- model.sel(mod_list, rank="DIC")

model_selection
(Intercept) inter_treatment mean.spl inter_treatment:mean.spl family df logLik DIC delta weight
mean.gaps ~ inter_treatment * mean.spl -0.0181886 + 0.0067237 + (NA) 7 69.40951 -127.7444 0.0000000 0.3578651599
mean.gaps ~ inter_treatment 0.6695200 + NA NA (NA) 5 68.40630 -127.6362 0.1081806 0.3390223244
mean.gaps ~ inter_treatment + mean.spl 0.2087143 + 0.0044737 NA (NA) 6 68.92134 -127.3879 0.3564624 0.2994432988
mean.gaps ~ 1 0.7031442 NA NA NA (NA) 4 63.05805 -118.0621 9.6822948 0.0028264143
mean.gaps ~ mean.spl 0.6270569 NA 0.0007432 NA (NA) 5 62.37306 -115.6420 12.1023576 0.0008428026

Best model summary

# summary
best_mod_gap <- gap_mods[[which(names(gap_mods) == row.names(model_selection)[1])]][[1]]

sm_gap <- as.data.frame(summary(best_mod_gap)$solutions[, -5])

sm_gap
post.mean l-95% CI u-95% CI eff.samp
(Intercept) -0.0181886 -0.7200648 0.7339871 2975.767
inter_treatmentBefore 0.5859397 -0.1818819 1.3476811 9700.000
mean.spl 0.0067237 -0.0000158 0.0140674 3457.288
inter_treatmentBefore:mean.spl -0.0051698 -0.0129531 0.0022992 9700.000
dat_gaps <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")

# subset data
dat_gaps$ID <- as.character(dat_gaps$ID)

ggplot(dat_gaps, aes(x = mean.gaps, y = xcrr)) +
  geom_point(size = 4, color = viridis(5, alpha = 0.5)[3]) +
  labs(y = "Mean gap duration (s)", x = "Sound pressure level (dB)") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24) 

Effect size posterior distribution

# stack posteriors
Y <- stack(as.data.frame(best_mod_gap$Sol[, -1]))

# plot posteriors
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution",x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 0.9, hjust=1)) 

# zoom in
Y <- stack(as.data.frame(best_mod_gap$Sol[, -c(1, 2)]))
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution",x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 0.9, hjust=1)) 

  • No detectable effect on gap duration

 

Diagnostic plots for MCMCglmm models

Include Gelman/Rubin’s convergence diagnostic, MCMC chain trace (all tree replicates in a single plot: yellow, blue and red colors) and autocorrelation plots:

SPL and period of the day

mcmc_diagnostics(period_mods)
## [1] "model: Null"

## [1] "model: period"

SPL and condition

mcmc_diagnostics(size_mods)
## [1] "model: Null"

## [1] "model: PC1"

## [1] "model: lift.weight"

## [1] "model: lift.weight + PC1"

## [1] "model: lift.weight *  PC1"

SPL and interactions

mcmc_diagnostics(inter_mods)
## [1] "model: Null"

## [1] "model: inter_treatment"

## [1] "model: inter_treatment + inter_type"

## [1] "model: inter_treatment * inter_type"

SPL and spectrographic consistency

mcmc_diagnostics(song_consistency_mods)
## [1] "model: xcrr ~ 1"

## [1] "model: xcrr ~ treatment"

## [1] "model: xcrr ~ treatment + mean.spl"

## [1] "model: xcrr ~ treatment * mean.spl"

## [1] "model: xcrr ~mean.spl"

SPL and gap duration

mcmc_diagnostics(gap_mods)
## [1] "model: mean.gaps ~ 1"

## [1] "model: mean.gaps ~ inter_treatment"

## [1] "model: mean.gaps ~ inter_treatment + mean.spl"

## [1] "model: mean.gaps ~ inter_treatment * mean.spl"

## [1] "model: mean.gaps ~ mean.spl"


Session information

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=pt_PT.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_PT.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_PT.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] gridExtra_2.3      lme4_1.1-26        corrplot_0.84      MuMIn_1.43.17     
##  [5] MCMCglmm_2.29      ape_5.4-1          coda_0.19-4        Matrix_1.2-18     
##  [9] rptR_0.9.22        readxl_1.3.1       viridis_0.5.1      viridisLite_0.3.0 
## [13] Rraven_1.0.10      warbleR_1.1.27     NatureSounds_1.0.3 knitr_1.31        
## [17] seewave_2.1.6      tuneR_1.3.3        ggplot2_3.3.2      pbapply_1.4-3     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6        lattice_0.20-41   corpcor_1.6.9     fftw_1.0-6       
##  [5] digest_0.6.27     R6_2.5.0          cellranger_1.1.0  signal_0.7-6     
##  [9] stats4_4.0.3      evaluate_0.14     highr_0.8         pillar_1.4.6     
## [13] rlang_0.4.10      cubature_2.0.4.1  minqa_1.2.4       nloptr_1.2.2.2   
## [17] rmarkdown_2.4     labeling_0.3      splines_4.0.3     statmod_1.4.35   
## [21] stringr_1.4.0     RCurl_1.98-1.3    munsell_0.5.0     proxy_0.4-25     
## [25] compiler_4.0.3    xfun_0.22         pkgconfig_2.0.3   htmltools_0.5.1.1
## [29] tidyselect_1.1.0  tibble_3.0.4      tensorA_0.36.1    dtw_1.22-3       
## [33] crayon_1.4.1      dplyr_1.0.2       withr_2.3.0       MASS_7.3-53      
## [37] bitops_1.0-6      nlme_3.1-149      gtable_0.3.0      lifecycle_0.2.0  
## [41] magrittr_2.0.1    scales_1.1.1      stringi_1.5.3     farver_2.0.3     
## [45] ellipsis_0.3.1    generics_0.1.0    vctrs_0.3.6       boot_1.3-25      
## [49] rjson_0.2.20      tools_4.0.3       glue_1.4.2        purrr_0.3.4      
## [53] yaml_2.2.1        colorspace_1.4-1